#include "Riostream.h"
#include <iostream>
#include <iomanip>
#include <string>

#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TMath.h"
#include "TRandom.h"
#include "TFormula.h"
#include "TSystem.h"
#include "TEnv.h"
#include "TF1.h"
#include "TF2.h"
#include "TF3.h"

#include "CommandLineInterface.hh"
#include "Barrel.hh"
#include "Annular.hh"
#include "Germanium.hh"
#include "SRIMloader.hh"

using namespace TMath;
using namespace std;

fstream Nfile;

ClassImp(Barrel);
ClassImp(Annular);
ClassImp(Germanium);

SRIM srim;

TF2* r_fun = new TF2("r_fun","sqrt(pow([0]+x,2)+pow([1]+y,2))",-10.,10,-10.,10.);
TF3* R_fun = new TF3("R_fun","sqrt(pow([0]+x,2)+pow([1]+y,2)+pow([2]+z,2))",-10.,10,-10.,10.,-5.,5.);


Double_t AverageX(Double_t *x, Double_t *par){
	const Double_t pi = TMath::Pi();
	Double_t x_Nuc, x_par, y_Nuc, y_par;
	Double_t phi_Nuc, phi_par;
	Double_t time[5][5];

	x_par = par[0];
	y_par = par[1];
	x_Nuc = par[2];
	y_Nuc = par[3];
	
	if (x_par>0 && y_par>0) phi_par = atan((y_par-x[1])/(x_par-x[0]));
	if (x_par>0 && y_par<0) phi_par = atan((y_par-x[1])/(x_par-x[0]))+2*pi;
	if (x_par<0 && y_par>0) phi_par = atan((y_par-x[1])/(x_par-x[0]))+pi;
	if (x_par<0 && y_par<0) phi_par = atan((y_par-x[1])/(x_par-x[0]))+pi;
	if (x_Nuc>0 && y_Nuc>0) phi_Nuc = atan((y_Nuc-x[1])/(x_Nuc-x[0]));
	if (x_Nuc>0 && y_Nuc<0) phi_Nuc = atan((y_Nuc-x[1])/(x_Nuc-x[0]))+2*pi;
	if (x_Nuc<0 && y_Nuc>0) phi_Nuc = atan((y_Nuc-x[1])/(x_Nuc-x[0]))+pi;
	if (x_Nuc<0 && y_Nuc<0) phi_Nuc = atan((y_Nuc-x[1])/(x_Nuc-x[0]))+pi;
	return 0.5*(cos(phi_par)+cos(phi_Nuc));
}

Double_t AverageY(Double_t *x, Double_t *par){
	const Double_t pi = TMath::Pi();
	Double_t x_Nuc, x_par, y_Nuc, y_par;
	Double_t phi_Nuc, phi_par;
	
	x_par = par[0];
	y_par = par[1];
	x_Nuc = par[2];
	y_Nuc = par[3];
	
	if (x_par>0 && y_par>0) phi_par = atan((y_par-x[1])/(x_par-x[0]));
	if (x_par>0 && y_par<0) phi_par = atan((y_par-x[1])/(x_par-x[0]))+2*pi;
	if (x_par<0 && y_par>0) phi_par = atan((y_par-x[1])/(x_par-x[0]))+pi;
	if (x_par<0 && y_par<0) phi_par = atan((y_par-x[1])/(x_par-x[0]))+pi;
	if (x_Nuc>0 && y_Nuc>0) phi_Nuc = atan((y_Nuc-x[1])/(x_Nuc-x[0]));
	if (x_Nuc>0 && y_Nuc<0) phi_Nuc = atan((y_Nuc-x[1])/(x_Nuc-x[0]))+2*pi;
	if (x_Nuc<0 && y_Nuc>0) phi_Nuc = atan((y_Nuc-x[1])/(x_Nuc-x[0]))+pi;
	if (x_Nuc<0 && y_Nuc<0) phi_Nuc = atan((y_Nuc-x[1])/(x_Nuc-x[0]))+pi;

	return 0.5*(sin(phi_par)+sin(phi_Nuc));
}

Double_t H(Double_t *x, Double_t *par){
	Double_t x_Nuc, x_par, y_Nuc, y_par;
	
	x_par = par[0];
	y_par = par[1];
	x_Nuc = par[2];
	y_Nuc = par[3];
	
	return ((y_Nuc-y_par)*x[0] + (y_Nuc*x_par-y_par*x_Nuc))/(x_Nuc-x_par);
}

TF1* h = new TF1("h",H,-10.,10.,4);

Double_t G(Double_t *var, Double_t *par){
	Double_t x_Nuc, x_par, y_Nuc, y_par, z_Nuc, z_par;
	Double_t P_Nuc, P_par, P_in;
	Double_t x, y, z;
	Double_t R_Nuc, R_par;
	Double_t g1, g2;

	x_par = par[0];
	y_par = par[1];
	z_par = par[2];
	P_par = par[3];
	x_Nuc = par[4];
	y_Nuc = par[5];
	z_Nuc = par[6];
	P_Nuc = par[7];
	P_in = par[8];

//	x = par[9];
	y = par[10];
//	z = var[0];
	z = 0;
	x = var[0];
	h->SetParameters(x_par,y_par,x_Nuc,y_Nuc);
	y = h->Eval(x);

	R_fun->SetParameters(x_Nuc,y_Nuc,z_Nuc);
	R_Nuc = R_fun->Eval(x,y,z);
	R_fun->SetParameters(x_par,y_par,z_par);
	R_par = R_fun->Eval(x,y,z);

	g1 = (x_par+x)*(x_Nuc+x) + (y_par+y)*(y_Nuc+y) + (z_par+z)*(z_Nuc+z);
	g2 = 0.5*R_Nuc*R_par*((P_in*P_in-P_Nuc*P_Nuc-P_par*P_par)/(P_Nuc*P_par));

	return TMath::Abs(g1-g2); 
}

TF1* g = new TF1("g",G,-5.,5.,11);

Double_t Fthetain(Double_t *var, Double_t *par){
	Double_t x_Nuc, x_par, y_Nuc, y_par, z_Nuc, z_par;
	Double_t P_Nuc, P_par, P_in;
	Double_t x, y, z;
	Double_t r_Nuc, r_par, R_Nuc, R_par;
	Double_t f1, f2;

	x_par = par[0];
	y_par = par[1];
	z_par = par[2];
	P_par = par[3];
	x_Nuc = par[4];
	y_Nuc = par[5];
	z_Nuc = par[6];
	P_Nuc = par[7];
	P_in = par[8];

	y = par[9];
	x = var[0];

//	h->SetParameters(x_par,y_par,x_Nuc,y_Nuc);
//	y = h->Eval(x);
//	g->SetParameters(x_par,y_par,z_par,P_par,x_Nuc,y_Nuc,z_Nuc,P_Nuc,P_in,x,y);
//	z = g->GetMinimumX(-5.,5.);
	z = 0;

	r_fun->SetParameters(x_Nuc,y_Nuc);
	r_Nuc = r_fun->Eval(x,y);
	r_fun->SetParameters(x_par,y_par);
	r_par = r_fun->Eval(x,y);
	R_fun->SetParameters(x_Nuc,y_Nuc,z_Nuc);
	R_Nuc = R_fun->Eval(x,y,z);
	R_fun->SetParameters(x_par,y_par,z_par);
	R_par = R_fun->Eval(x,y,z);
	
	f1 = pow(P_par*r_par/R_par,2) + pow(P_Nuc*r_Nuc/R_Nuc,2) + 2*P_par*P_Nuc/(R_par*R_Nuc)*((x_par+x)*(x_Nuc+x)+(y_par+y)*(y_Nuc+y));
	f2 = P_par*(z_par+z)/R_par + P_Nuc*(z_Nuc+z)/R_Nuc;

	return sqrt(f1)/f2;

}

TF1* fthetain = new TF1("fthetain",Fthetain,-10.,10.,9);


int main(int argc, char* argv[]){
	vector<char*> InputFiles;
	vector<char*> Histos;
	char* OutputFile = NULL;
	char* MBFile = NULL;
	char* cutfile = NULL;
	char* srimdir = NULL;
	char* SetFile = NULL;
	bool verbose = false;
	int tiltcorr = 0;
	int sim = 0;
	
	const Double_t pi = TMath::Pi();

	Double_t m[150];
	Double_t Ein, Ein_set, theta_in, phi_in, xin, yin, zin;
	Double_t P_in, P_in_x, P_in_y, P_in_z;
	Double_t Pxin, Pyin, Pzin;
	Double_t ttarget, Rafter;
	Double_t Q, Eex;
	Double_t CDThreshold[4];
	Int_t Apart, Atarget;
	Double_t xpad[4], ypad[4], zpad[4];
	
	Double_t PT, Edet[3], dEcor, Ecor;
	Double_t E_par, E_Nuc, P_par, P_Nuc, theta_par;
	Double_t P_par_x, P_par_y, P_par_z, P_Nuc_x, P_Nuc_y, P_Nuc_z;
	Double_t R_par, R_Nuc;
	Double_t time[5][5];

	Double_t x, y, z, r, alpha;
	Double_t theta, phi, thetadet;
	Double_t x_par, y_par, z_par, x_Nuc, y_Nuc, z_Nuc;
	Double_t x1, x2;

	// Masses
	m[1] = 1.00782503207*931.494e3;
	m[2] = 2.01410177785*931.494e3;
	m[3] = 3.01604927767*931.494e3;
	m[4] = 4.00260325415*931.494e3;
	m[9] = 9.012182201*931.494e3;
	m[10] = 10.013533818*931.494e3;
	m[11] = 11.021657749*931.494e3;
	m[12] = 12.026920737*931.494e3; // 12Be
	m[20] = 19.999981315*931.494e3; // 20F
	m[21] = 20.999948951*931.494e3; // 21F
	m[22] = 22.002998815*931.494e3; // 22F
	m[23] = 22.994466904*931.494e3; // 23Ne
	m[107] = 106.905096820*931.494e3;

	CommandLineInterface* interface = new CommandLineInterface();

	interface->Add("-i", "inputfiles", &InputFiles);
	interface->Add("-o", "outputfile", &OutputFile);
	interface->Add("-M", "miniball angle", &MBFile);
	interface->Add("-C", "cut functions", &cutfile);
	interface->Add("-S", "SRIM files directory", &srimdir);
	interface->Add("-Set", "Setting file", &SetFile);
	interface->Add("-t", "1 tiltcorr histos, 0 no tiltcorr", &tiltcorr);
	interface->Add("-v", "verbose", &verbose);
	interface->Add("-sim", "Set if you are analyzing a simulation",&sim);


	interface->CheckFlags(argc, argv);

	if(InputFiles.size() == 0 || OutputFile == NULL){
		cerr<<"You have to provide at least one input file and the output file!"<<endl;
		exit(1);
	}
	if(srimdir == NULL){
		cerr << "SRIM directory is needed!"<<endl;
		exit(1);
	}

	cout<<"input file(s):"<<endl;
	for(unsigned int i=0; i<InputFiles.size(); i++){
		cout<<InputFiles[i]<<endl;
	}

	cout<<"output file: "<<OutputFile<< endl;

	// Loading settings
	TEnv* Settings = new TEnv(SetFile);
	
	// Beam energy
	Ein_set = Settings->GetValue("E.beam",2.85e3);
	Ein_set = Ein_set*11;
	theta_in = Settings->GetValue("theta.beam",0.);
	theta_in = theta_in/180.*pi;
	phi_in = Settings->GetValue("phi.beam",0.);
	phi_in = phi_in/180.*pi;
	xin = Settings->GetValue("x.beam",0);
	yin = Settings->GetValue("y.beam",0.);

	ttarget = Settings->GetValue("t.target",10.);
	Atarget = Settings->GetValue("A.target",2);

	if (Atarget==12) m[12] = 12.000000000*931.494e3; // 12C

	// Pad settings
	for (Int_t ipad=0;ipad<4;ipad++){
		xpad[ipad] = Settings->GetValue(Form("x.pad.%d",ipad),0);
		ypad[ipad] = Settings->GetValue(Form("y.pad.%d",ipad),0);
		zpad[ipad] = Settings->GetValue(Form("z.pad.%d",ipad),0);
	}

	// A for particles stopped in dE
	Apart = Settings->GetValue("A.particle",2);

	// Generating SRIM graphs for energyloss calculations
	srim.Creategraphs(srimdir);
	//printf("R = %f, E = %f \n",srim.RMylar[4]->Eval(3.0),srim.EMylar[4]->Eval(11.0));
	if (sim==0){
	Rafter = srim.RT[11]->Eval(Ein_set);
	Ein_set = srim.ET[11]->Eval(Rafter - ttarget/(2.*cos(theta_in)));
	}
	
	Pxin = sqrt(2*m[11]*Ein_set)*sin(theta_in)*cos(phi_in);
	Pyin = sqrt(2*m[11]*Ein_set)*sin(theta_in)*sin(phi_in);
	Pzin = sqrt(2*m[11]*Ein_set)*cos(theta_in);
	
	


	TChain* tr;
	tr = new TChain("caltr");
	for(unsigned int i=0; i<InputFiles.size(); i++){
		tr->Add(InputFiles[i]);
	}

	if(tr == NULL){
		cout << "could not find tree caltr in file " << endl;
		for(unsigned int i=0; i<InputFiles.size(); i++){
			cout<<InputFiles[i]<<endl;
		}
		return 3;
	}

	vector<Barrel> *FBarrel = new vector<Barrel>;
	vector<Barrel> *BBarrel = new vector<Barrel>;
	vector<Annular> *FCD = new vector<Annular>;
	vector<Germanium> *Miniball = new vector<Germanium>;
	long long EbisTime;
	long long T1Time;
	long long SuperCycleTime;

	tr->SetBranchAddress("ForwardBarrel",&FBarrel);
	tr->SetBranchAddress("BackwardBarrel",&BBarrel);
	tr->SetBranchAddress("ForwardCD",&FCD);
	tr->SetBranchAddress("Miniball",&Miniball);
	tr->SetBranchAddress("EbisTime",&EbisTime);
	tr->SetBranchAddress("T1Time",&T1Time);
	tr->SetBranchAddress("SuperCycleTime",&SuperCycleTime);

	TFile* outfile = new TFile(OutputFile,"recreate");

	if(outfile->IsZombie()){
		return 4;
	}

	CDThreshold[0] = 500.;
	CDThreshold[1] = 1300.;
	CDThreshold[2] = 500.;
	CDThreshold[3] = 800.;

	// Generating cut functions
	TFormula* For_cut[4][16];
	TFormula* Back_cut[4];
	TFormula* CD_cut[6][4];
	TFormula* CD_PT_cut[6][4];
	TFormula* Be_cut[4];

	TEnv *Cutfunc = new TEnv(cutfile);

	// Forward Barrel
	for (Int_t istrip=0;istrip<16;istrip++){
		for (Int_t icut=0;icut<4;icut++) For_cut[icut][istrip] = new TFormula(Form("For_cut_%d_%d",icut,istrip),Cutfunc->GetValue(Form("For.%d.%d",icut,istrip),"0*x"));
	}
	// 11 Beryllium
	Be_cut[0] = new TFormula("Be_cut_0",Cutfunc->GetValue("For.11","0*x"));
	Be_cut[1] = new TFormula("Be_cut_1",Cutfunc->GetValue("For.12","0*x"));

	// Backward Barrel
	for (Int_t icut=0;icut<4;icut++) Back_cut[icut] = new TFormula(Form("Back_cut_%d",icut),Cutfunc->GetValue(Form("Back.%d",icut),"0*x"));

	// CD
	for (Int_t idet=0;idet<4;idet++){
		for (Int_t icut=0;icut<6;icut++){
			CD_cut[icut][idet] = new TFormula(Form("CD_cut_%d_%d",icut,idet),Cutfunc->GetValue(Form("CD.%d.%d",icut,idet),"0*x"));
			CD_PT_cut[icut][idet] = new TFormula(Form("CD_PT_cut_%d_%d",icut,idet),Cutfunc->GetValue(Form("CD_PT.%d.%d",icut,idet),"0*x"));
		}
	}
	// 11 Beryllium
	Be_cut[2] = new TFormula("Be_cut_2",Cutfunc->GetValue("CD.11","0*x"));
	Be_cut[3] = new TFormula("Be_cut_3",Cutfunc->GetValue("CD.12","0*x"));

	
	TF2* fAveragePhix = new TF2("fAveragePhix",AverageX,-20,20,-20,20,4);
	TF2* fAveragePhiy = new TF2("fAveragePhiy",AverageY,-20,20,-20,20,4);


	// Prepare for vectorial histograms
	TH2F* Triton_XY_det[4];
	TH2F* Triton_Thetain_det[4];
	TH2F* Triton_XvsThetain_det[4];
	TH2F* Triton_YvsThetain_det[4];
	TH2F* Triton_XvsdeltaEin1_det[4];
	TH2F* Triton_YvsdeltaEin1_det[4];
	TH2F* Triton_XvsdeltaEin2_det[4];
	TH2F* Triton_YvsdeltaEin2_det[4];

	// Deuterons
	TH2F* Deuteron_Average = new TH2F("Deuteron_Average","Phi averag for t+10Be",200,-1,1,200,-1,1);
	TH2F* Deuteron_Average_XY = new TH2F("Deuteron_Average_XY","Phi average with new orrigo for t+10Be",200,-1,1,200,-1,1);
	TH2F* Deuteron_XY = new TH2F("Deuteron_XY","X vs Y for t+10Be",200,-10,10,200,-10,10);
	TH1F* Deuteron_Ein1 = new TH1F("Deuteron_Ein1","Ein for t+10Be",400,0,40);
	TH1F* Deuteron_Ein2 = new TH1F("Deuteron_Ein2","Ein for t+10Be",400,0,40);
	
	TH2F* Deuteron_dEvsEcor = new TH2F("Deuteron_dEvsEcor","dE vs E",200,0,30,200,0,15);
	TH2F* Deuteron_Thetain = new TH2F("Deuteron_Thetain","Theta in with cos(phi) vs sin(phi)",200,-200,200,200,-200,200);
	TH2F* Deuteron_XvsThetain = new TH2F("Deuteron_XvsThetain","x vs Theta*cos(phi)",200,-10,10,200,-200,200);
	TH2F* Deuteron_YvsThetain = new TH2F("Deuteron_YvsThetain","y vs Theta*sin(phi)",200,-10,10,200,-200,200);
	TH2F* Deuteron_EvsTheta = new TH2F("Deuteron_EvsTheta","theta vs E for t",180,0,90,350,0,35);
	TH2F* Be11_EvsTheta = new TH2F("Be11_EvsTheta","theta vs E for 10Be",180,0,90,350,0,35);
	
	// Tritons
	TH2F* Triton_Average = new TH2F("Triton_Average","Phi averag for t+10Be",200,-1,1,200,-1,1);
	TH2F* Triton_Average_u2 = new TH2F("Triton_Average_u2","Phi average for t+10Be",200,-.1,.1,200,-.1,.1);
	TH2F* Triton_Average_XY = new TH2F("Triton_Average_XY","Phi average with new orrigo for t+10Be",200,-1,1,200,-1,1);
	TH2F* Triton_Average_XY_2 = new TH2F("Triton_Average_XY_2","Phi average with new orrigo for t+10Be",200,-1,1,200,-1,1);
	TH2F* Triton_XY = new TH2F("Triton_XY","X vs Y for t+10Be",200,-10,10,200,-10,10);
	TH2F* Triton_XY_2 = new TH2F("Triton_XY_2","X vs Y for t+10Be",200,-10,10,200,-10,10);
	TH2F* Triton_XY_3 = new TH2F("Triton_XY_3","X vs Y for t+10Be",200,-10,10,200,-10,10);
	TH2F* Triton_XY_Strip = new TH2F("Triton_XY_Strip","X vs Y for t+10Be",200,-10,10,200,-10,10);
	TH1F* Triton_Z = new TH1F("Triton_Z","Z for t+10Be",200,-5,5);
	TH1F* Triton_Ein1 = new TH1F("Triton_Ein1","Ein for t+10Be",400,0,40);
	TH1F* Triton_Ein2 = new TH1F("Triton_Ein2","Ein for t+10Be",400,0,40);
	TH1F* Triton_Eex = new TH1F("Triton_Eex","E* for 10Be",1000,-2,8.);
	TH2F* Triton_ZXpar = new TH2F("Triton_ZXpar","X_par vs Z_par for t+10Be",200,0,60,200,-26,26);
	TH2F* Triton_ZYpar = new TH2F("Triton_ZYpar","X_par vs Z_par for t+10Be",200,0,60,200,-26,26);
	
	TH2F* For_dEvsEcor = new TH2F("For_dEvsEcor","dE vs E",300,0,30,300,0,15);
	TH2F* Triton_dEvsEcor = new TH2F("Triton_dEvsEcor","dE vs E",200,0,30,200,0,15);
	TH2F* Triton_Thetain = new TH2F("Triton_Thetain","Theta in with cos(phi) vs sin(phi)",200,-200,200,200,-200,200);
	TH2F* Triton_Thetain_XY = new TH2F("Triton_Thetain_XY","Theta in with cos(phi) vs sin(phi)",200,-20,20,200,-20,20);
	TH1F* Triton_Thetain_2 = new TH1F("Triton_Thetain_2","Theta",200,-200,200);
	TH2F* Triton_Thetain_Phiin_2 = new TH2F("Triton_Thetain_Phiin_2","Theta in with cos(phi) vs sin(phi)",200,-200,200,200,-200,200);
	TH2F* Triton_XvsThetain = new TH2F("Triton_XvsThetain","x vs Theta*cos(phi)",200,-10,10,200,-200,200);
	TH2F* Triton_YvsThetain = new TH2F("Triton_YvsThetain","y vs Theta*sin(phi)",200,-10,10,200,-200,200);
	TH2F* Triton_XvsdeltaEin1 = new TH2F("Triton_XvsdeltaEin1","x vs Ein-Ein_set",200,-10,10,200,-10,10);
	TH2F* Triton_YvsdeltaEin1 = new TH2F("Triton_YvsdeltaEin1","y vs Ein-Ein_set",200,-10,10,200,-10,10);
	TH2F* Triton_XvsdeltaEin2 = new TH2F("Triton_XvsdeltaEin2","x vs Ein-Ein_set",200,-10,10,200,-10,10);
	TH2F* Triton_YvsdeltaEin2 = new TH2F("Triton_YvsdeltaEin2","y vs Ein-Ein_set",200,-10,10,200,-10,10);
	
	TH2F* Triton_dtvsEex = new TH2F("Triton_dtvsEex","hmm",100,-1250,1250,500,-2,8);

	TH2F* Triton_EvsTheta = new TH2F("Triton_EvsTheta","theta vs E for t",180,0,90,350,0,35);
	TH2F* Be10_EvsTheta = new TH2F("Be10_EvsTheta","theta vs E for 10Be",180,0,90,350,0,35);

	for (Int_t idet=0;idet<4;idet++){
		Triton_XY_det[idet] = new TH2F(Form("Triton_XY_det_%d",idet),Form("X vs Y in det %d",idet),200,-10,10,200,-10,10);
		Triton_Thetain_det[idet] = new TH2F(Form("Triton_Thetain_det_%d",idet),"Theta in with cos(phi) vs sin(phi)",200,-200,200,200,-200,200);
		Triton_XvsThetain_det[idet] = new TH2F(Form("Triton_XvsThetain_det_%d",idet),"y vs Theta*sin(phi)",200,-10,10,200,-200,200);
		Triton_YvsThetain_det[idet] = new TH2F(Form("Triton_YvsThetain_det_%d",idet),"y vs Theta*sin(phi)",200,-10,10,200,-200,200);
		Triton_XvsdeltaEin1_det[idet] = new TH2F(Form("Triton_XvsdeltaEin1_det_%d",idet),"x vs Ein-Ein_set",200,-10,10,200,-10,10);
		Triton_YvsdeltaEin1_det[idet] = new TH2F(Form("Triton_YvsdeltaEin1_det_%d",idet),"y vs Ein-Ein_set",200,-10,10,200,-10,10);
		Triton_XvsdeltaEin2_det[idet] = new TH2F(Form("Triton_XvsdeltaEin2_det_%d",idet),"x vs Ein-Ein_set",200,-10,10,200,-10,10);
		Triton_YvsdeltaEin2_det[idet] = new TH2F(Form("Triton_YvsdeltaEin2_det_%d",idet),"y vs Ein-Ein_set",200,-10,10,200,-10,10);
	}
	
	// END HISTOGRAMS

	bool partfound;
	Double_t nentries = tr->GetEntries();
	Int_t nbytes = 0;
	Int_t status;

	for(int i=0; i<nentries;i++){
		
		status = tr->GetEvent(i);
		if(status == -1){
			cerr<<"Error occured, couldn't read entry "<<i<<" from tree "<<tr->GetName()<<" in file "<<tr->GetFile()->GetName()<<endl;
			return 5;
		}
		else if(status == 0){
			cerr<<"Error occured, entry "<<i<<" in tree "<<tr->GetName()<<" in file "<<tr->GetFile()->GetName()<<" doesn't exist"<<endl;
			return 6;
		}
		nbytes += status;


		partfound =false;

		// Forward barrel
		for(unsigned int j_bar=0; j_bar<FBarrel->size(); j_bar++){
			if((*FBarrel)[j_bar].GetMult() ==0) continue;

			if( (*FBarrel)[j_bar].GetID() < 0 || (*FBarrel)[j_bar].GetID() > 3 ){
				cerr<<"Error in entry "<<i<<": "<<j_bar<<". forward barrel detector id is wrong: "<< (*FBarrel)[j_bar].GetID()<<endl;
				continue;
			}

			if( (*FBarrel)[j_bar].GetStripNr()[0] < 0 || (*FBarrel)[j_bar].GetStripNr()[0] > 15 ){
				cerr<<"Error in entry "<<i<<": "<<(*FBarrel)[j_bar].GetID()<<". forward barrel detector strip nr is wrong: "<<(*FBarrel)[j_bar].GetStripNr()[0]<<endl;
				continue;
			}
			
			partfound=true;
			time[0][j_bar] = (*FBarrel)[j_bar].GetTime();
			Edet[0] = (*FBarrel)[j_bar].GetEdet();
			
			// Cylindrical coordinates
			x = -50*((*FBarrel)[j_bar].GetStripPos()[0]-.5) + xin*cos(pi/2.*(*FBarrel)[j_bar].GetID()) + yin*sin(pi/2.*(*FBarrel)[j_bar].GetID());
			y = 29. + xin*sin(pi/2.*(*FBarrel)[j_bar].GetID()) - yin*cos(pi/2.*(*FBarrel)[j_bar].GetID());
			z = 8. + 3.125*((*FBarrel)[j_bar].GetStripNr()[0] + gRandom->Rndm()); 
			// Moving each detector
			for (Int_t ipad=0;ipad<4;ipad++){
				if ((*FBarrel)[j_bar].GetID()==ipad){
					x += xpad[ipad];
					y += ypad[ipad];
					z += zpad[ipad];
				}
			}

			r = sqrt(pow(y,2) + pow(x,2));
			alpha = atan(x/y);
			// Spherical coordinates
			theta = atan(r/z);
			phi = alpha + pi/2.*(-(*FBarrel)[j_bar].GetID()+1);
			if (phi< 0) phi += 2*pi;
			if (phi> 2*pi) phi -= 2*pi;
			// Detector angle
			thetadet = acos(sin(theta)*cos(alpha));
			
			// Filling histos
			
			
			// No punchthrough
			PT = srim.EDL[Apart]->Eval(140./cos(thetadet));
			PT = 0.;


			// Particle identification
			if (((*FBarrel)[j_bar].GetEdet() > 400.) && ((*FBarrel)[j_bar].GetRear() > PT)){
				Ecor = (*FBarrel)[j_bar].GetEdet() + (1. - cos(thetadet))*(*FBarrel)[j_bar].GetRear();
				dEcor = (*FBarrel)[j_bar].GetRear()*cos(thetadet);
				

				if ((*FBarrel)[j_bar].GetID()!=1) For_dEvsEcor->Fill(Ecor/1.e3,dEcor/1.e3);

				// Deuterons
				if ((dEcor > For_cut[1][(*FBarrel)[j_bar].GetStripNr()[0]]->Eval(Ecor)) && (dEcor < For_cut[2][(*FBarrel)[j_bar].GetStripNr()[0]]->Eval(Ecor))){
					// Position
					x_par = r*cos(phi);
					y_par = r*sin(phi);
					z_par = z;

					Deuteron_dEvsEcor->Fill(Ecor/1.e3,dEcor/1.e3);

					// Energy and momentum
					Rafter = srim.RMylar[2]->Eval((*FBarrel)[j_bar].GetEdet() + (*FBarrel)[j_bar].GetRear());
					E_par = srim.EMylar[2]->Eval(Rafter + 11.57/cos(thetadet));
					Rafter = srim.RT[2]->Eval(E_par);
					E_par = srim.ET[2]->Eval(Rafter + ttarget/(2.*cos(theta)));
					
//					Rafter = srim.RT[2]->Eval(E_par);
//					E_par = srim.ET[2]->Eval(Rafter + ttarget/(2.*cos(theta)));
					P_par = sqrt(2*m[2]*E_par);	
					P_par_x = P_par*sin(theta)*cos(phi); 
					P_par_y = P_par*sin(theta)*sin(phi); 
					P_par_z = P_par*cos(theta); 
					theta_par = theta;
					
					
					// Forward CD
					for(unsigned int j_CD=0; j_CD<FCD->size(); j_CD++){
						if( (*FCD)[j_CD].GetID() < 0 || (*FCD)[j_CD].GetID() >= 4 ){
							cerr<<"Error in entry "<<i<<": "<<j_CD<<". backward barrel detector id is wrong: "<<(*FCD)[j_CD].GetID()<<endl;
							continue;
						}
					
						// Ring run
						for( unsigned int k=0;k<(*FCD)[j_CD].GetRingNr().size();k++){
							// Strip run
							for( unsigned int l=0;l<(*FCD)[j_CD].GetStripNr().size();l++){
								// Real events
								if (TMath::Abs((*FCD)[j_CD].GetRingEnergy()[k] - (*FCD)[j_CD].GetStripEnergy()[l]) < 500){
									
									// Angle
									// Values from centrum
									r = 9. + (gRandom->Rndm() + (*FCD)[j_CD].GetRingNr()[k])*2.;				
									z_Nuc = 63.;
									theta = atan(r/z_Nuc);
						
									if (((*FCD)[j_CD].GetStripNr()[l]>=4)&&((*FCD)[j_CD].GetStripNr()[l]<12)){
										phi = -.0593411*(4+ 2*((*FCD)[j_CD].GetStripNr()[l]-4 + gRandom->Rndm())) + .71209;
									}
									else if ((*FCD)[j_CD].GetStripNr()[l]<4){
										phi = -.0593411*((*FCD)[j_CD].GetStripNr()[l] + gRandom->Rndm()) + .71209;
									}
									else{
										phi = -.0593411*((*FCD)[j_CD].GetStripNr()[l]+8 + gRandom->Rndm()) + .71209;
									}
							
									phi = phi + pi/2.*(-(*FCD)[j_CD].GetID()+1);


									// Determine x and y coordinates
									x_Nuc = r*cos(phi);
									y_Nuc = r*sin(phi);							
								
									// Determining the energy
									Rafter = srim.RDL[11]->Eval((*FCD)[j_CD].GetRingEnergy()[k]);
									E_Nuc = srim.EDL[11]->Eval(Rafter + .7/cos(theta));
									Rafter = srim.RT[11]->Eval(E_Nuc);
									E_Nuc = srim.ET[11]->Eval(Rafter + ttarget/(2.*cos(theta)));
									P_Nuc = sqrt(2*m[11]*E_Nuc);
									P_Nuc_x = P_Nuc*sin(theta)*cos(phi); 
									P_Nuc_y = P_Nuc*sin(theta)*sin(phi); 
									P_Nuc_z = P_Nuc*cos(theta); 
									
									
									Deuteron_EvsTheta->Fill(theta_par/pi*180.,E_par/1.e3);
									Be11_EvsTheta->Fill(theta/pi*180.,E_Nuc/1.e3);
									
									// SIMPLE
									// Calculating the offset of the beam
//									if (((*FCD)[j_CD].GetStripNr()[l] < 4) || ((*FCD)[j_CD].GetStripNr()[l] > 12)) {
									
										fAveragePhix->SetParameters(x_par,y_par,x_Nuc,y_Nuc);
										fAveragePhiy->SetParameters(x_par,y_par,x_Nuc,y_Nuc);
										
										
//										x = ((x_par*y_Nuc - x_Nuc*y_par)*(y_par-y_Nuc))/((x_par-x_Nuc)*(x_Nuc-x_par)+(y_par-y_Nuc)*(y_Nuc-y_par));
//										y = ((x_par*y_Nuc - x_Nuc*y_par)*(x_Nuc-x_par))/((x_par-x_Nuc)*(x_Nuc-x_par)+(y_par-y_Nuc)*(y_Nuc-y_par));
										x = ((x_par*y_Nuc - x_Nuc*y_par)*(y_Nuc-y_par))/(pow((x_par-x_Nuc),2)+pow((y_par-y_Nuc),2));
										y = ((x_par*y_Nuc - x_Nuc*y_par)*(x_par-x_Nuc))/(pow((x_par-x_Nuc),2)+pow((y_par-y_Nuc),2));

										Deuteron_XY->Fill(x,y);
										Deuteron_Average->Fill(fAveragePhix->Eval(0,0),fAveragePhiy->Eval(0,0));
										Deuteron_Average_XY->Fill(fAveragePhix->Eval(x,y),fAveragePhiy->Eval(x,y));
//									}

									// Calculating the incoming momentum
									P_in_x = P_par_x+P_Nuc_x;	
									P_in_y = P_par_y+P_Nuc_y;	
									P_in_z = P_par_z+P_Nuc_z;	
								
									
									
									// Calculating the incoming angle
									theta_in = atan(sqrt(pow(P_in_x,2)+pow(P_in_y,2))/P_in_z);
									if (P_in_x>=0&&P_in_y>=0) phi_in = atan(P_in_y/P_in_x);
									if (P_in_x>=0&&P_in_y<0) phi_in = 2*pi + atan(P_in_y/P_in_x);
									if (P_in_x<0&&P_in_y>=0) phi_in = pi + atan(P_in_y/P_in_x);
									if (P_in_x<0&&P_in_y<0) phi_in = pi + atan(P_in_y/P_in_x);
									
									Deuteron_Thetain->Fill(theta_in*cos(phi_in)*1e3,theta_in*sin(phi_in)*1e3);
								
									Deuteron_XvsThetain->Fill(x,theta_in*cos(phi_in)*1e3);
									Deuteron_YvsThetain->Fill(y,theta_in*sin(phi_in)*1e3);
									Ein = (pow(P_in_x,2) + pow(P_in_y,2) + pow(P_in_z,2))/(2*m[11]);
									Deuteron_Ein2->Fill(Ein/1.e3);
									
									
									
									// ADVANCED
									// Calculating the offset
									
									Eex = 0;
									Q = m[11]+m[2]-m[11]-m[2];
									Ein = E_Nuc + E_par - Q - Eex;
									P_in = sqrt(2*m[11]*Ein);
									
									Deuteron_Ein1->Fill(Ein/1.e3);

								}
							}
						}
					}
				}	
				
				// Tritons
				if ((dEcor > For_cut[2][(*FBarrel)[j_bar].GetStripNr()[0]]->Eval(Ecor)) && (dEcor < For_cut[3][(*FBarrel)[j_bar].GetStripNr()[0]]->Eval(Ecor))){
					// Position
					x_par = r*cos(phi);
					y_par = r*sin(phi);
					z_par = z;

					Triton_dEvsEcor->Fill(Ecor/1.e3,dEcor/1.e3);

					// Energy and momentum
					Rafter = srim.RMylar[3]->Eval((*FBarrel)[j_bar].GetEdet() + (*FBarrel)[j_bar].GetRear());
					E_par = srim.EMylar[3]->Eval(Rafter + 11.57/cos(thetadet));
					Rafter = srim.RT[3]->Eval(E_par);
					E_par = srim.ET[3]->Eval(Rafter + ttarget/(2.*cos(theta)));
					
					P_par = sqrt(2*m[3]*E_par);	
					P_par_x = P_par*sin(theta)*cos(phi); 
					P_par_y = P_par*sin(theta)*sin(phi); 
					P_par_z = P_par*cos(theta); 
					theta_par = theta;
					
					
					// Forward CD
					for(unsigned int j_CD=0; j_CD<FCD->size(); j_CD++){
						if( (*FCD)[j_CD].GetID() < 0 || (*FCD)[j_CD].GetID() >= 4 ){
							cerr<<"Error in entry "<<i<<": "<<j_CD<<". backward barrel detector id is wrong: "<<(*FCD)[j_CD].GetID()<<endl;
							continue;
						}
					
						// Ring run
						for( unsigned int k=0;k<(*FCD)[j_CD].GetRingNr().size();k++){
							// Strip run
							for( unsigned int l=0;l<(*FCD)[j_CD].GetStripNr().size();l++){
								// Real events
								if (TMath::Abs((*FCD)[j_CD].GetRingEnergy()[k] - (*FCD)[j_CD].GetStripEnergy()[l]) < 500){
									time[1][j_CD] = (*FCD)[j_CD].GetTime();
									
									// Angle
									// Values from centrum
									r = 9. + (gRandom->Rndm() + (*FCD)[j_CD].GetRingNr()[k])*2.;				
									z_Nuc = 63.;
									theta = atan(r/z_Nuc);
						
									if (((*FCD)[j_CD].GetStripNr()[l]>=4)&&((*FCD)[j_CD].GetStripNr()[l]<12)){
										phi = -.0593411*(4+ 2*((*FCD)[j_CD].GetStripNr()[l]-4 + gRandom->Rndm())) + .71209;
									}
									else if ((*FCD)[j_CD].GetStripNr()[l]<4){
										phi = -.0593411*((*FCD)[j_CD].GetStripNr()[l] + gRandom->Rndm()) + .71209;
									}
									else{
										phi = -.0593411*((*FCD)[j_CD].GetStripNr()[l]+8 + gRandom->Rndm()) + .71209;
									}
							
									phi = phi + pi/2.*(-(*FCD)[j_CD].GetID()+1);


									// Determine x and y coordinates
									x_Nuc = r*cos(phi);
									y_Nuc = r*sin(phi);							
								
									// Determining the energy
									Rafter = srim.RDL[10]->Eval((*FCD)[j_CD].GetRingEnergy()[k]);
									E_Nuc = srim.EDL[10]->Eval(Rafter + .7/cos(theta));
									Rafter = srim.RT[10]->Eval(E_Nuc);
									E_Nuc = srim.ET[10]->Eval(Rafter + ttarget/(2.*cos(theta)));
									P_Nuc = sqrt(2*m[10]*E_Nuc);
									P_Nuc_x = P_Nuc*sin(theta)*cos(phi); 
									P_Nuc_y = P_Nuc*sin(theta)*sin(phi); 
									P_Nuc_z = P_Nuc*cos(theta); 
									
									Triton_dtvsEex->Fill((time[0][j_bar]-time[2][j_CD])*25,Eex/1000.);
									
									Triton_EvsTheta->Fill(theta_par/pi*180.,E_par/1.e3);
									Be10_EvsTheta->Fill(theta/pi*180.,E_Nuc/1.e3);
									
									// SIMPLE
									// Calculating the offset of the beam
//									if (((*FCD)[j_CD].GetStripNr()[l] < 4) || ((*FCD)[j_CD].GetStripNr()[l] > 12)) {
									
										fAveragePhix->SetParameters(x_par,y_par,x_Nuc,y_Nuc);
										fAveragePhiy->SetParameters(x_par,y_par,x_Nuc,y_Nuc);
										
										
//										x = ((x_par*y_Nuc - x_Nuc*y_par)*(y_par-y_Nuc))/((x_par-x_Nuc)*(x_Nuc-x_par)+(y_par-y_Nuc)*(y_Nuc-y_par));
//										y = ((x_par*y_Nuc - x_Nuc*y_par)*(x_Nuc-x_par))/((x_par-x_Nuc)*(x_Nuc-x_par)+(y_par-y_Nuc)*(y_Nuc-y_par));
										x = ((x_par*y_Nuc - x_Nuc*y_par)*(y_Nuc-y_par))/(pow((x_par-x_Nuc),2)+pow((y_par-y_Nuc),2));
										y = ((x_par*y_Nuc - x_Nuc*y_par)*(x_par-x_Nuc))/(pow((x_par-x_Nuc),2)+pow((y_par-y_Nuc),2));

										Triton_XY->Fill(x,y);
										Triton_XY_det[(*FBarrel)[j_bar].GetID()]->Fill(x,y);
										if (((*FCD)[j_CD].GetStripNr()[l]==8)&&((*FCD)[j_CD].GetID()==2)) Triton_XY_Strip->Fill(x,y);
										
										Triton_ZXpar->Fill(z_par,x_par);
										Triton_ZYpar->Fill(z_par,y_par);
										Triton_Average->Fill(fAveragePhix->Eval(0,0),fAveragePhiy->Eval(0,0));
										Triton_Average_u2->Fill(fAveragePhix->Eval(0,0),fAveragePhiy->Eval(0,0));
										Triton_Average_XY->Fill(fAveragePhix->Eval(x,y),fAveragePhiy->Eval(x,y));
//									}

									// Calculating the incoming momentum
									P_in_x = P_par_x+P_Nuc_x;	
									P_in_y = P_par_y+P_Nuc_y;	
									P_in_z = P_par_z+P_Nuc_z;	
								
									
									
									// Calculating the incoming angle
									theta_in = atan(sqrt(pow(P_in_x,2)+pow(P_in_y,2))/P_in_z);
									if (P_in_x>=0&&P_in_y>=0) phi_in = atan(P_in_y/P_in_x);
									if (P_in_x>=0&&P_in_y<0) phi_in = 2*pi + atan(P_in_y/P_in_x);
									if (P_in_x<0&&P_in_y>=0) phi_in = pi + atan(P_in_y/P_in_x);
									if (P_in_x<0&&P_in_y<0) phi_in = pi + atan(P_in_y/P_in_x);
									
									Triton_Thetain->Fill(theta_in*cos(phi_in)*1e3,theta_in*sin(phi_in)*1e3);
									Triton_Thetain_det[(*FBarrel)[j_bar].GetID()]->Fill(theta_in*cos(phi_in)*1e3,theta_in*sin(phi_in)*1e3);
								
									Triton_XvsThetain->Fill(x,theta_in*cos(phi_in)*1e3);
									Triton_YvsThetain->Fill(y,theta_in*sin(phi_in)*1e3);
									Triton_XvsThetain_det[(*FBarrel)[j_bar].GetID()]->Fill(x,theta_in*cos(phi_in)*1e3);
									Triton_YvsThetain_det[(*FBarrel)[j_bar].GetID()]->Fill(y,theta_in*sin(phi_in)*1e3);
									Ein = (pow(P_in_x,2) + pow(P_in_y,2) + pow(P_in_z,2))/(2*m[11]);
									Triton_Ein2->Fill(Ein/1.e3);
								
									Triton_XvsdeltaEin2->Fill((Ein-Ein_set)/1e3,x);
									Triton_YvsdeltaEin2->Fill((Ein-Ein_set)/1e3,y);
									Triton_XvsdeltaEin2_det[(*FBarrel)[j_bar].GetID()]->Fill((Ein-Ein_set)/1e3,x);
									Triton_YvsdeltaEin2_det[(*FBarrel)[j_bar].GetID()]->Fill((Ein-Ein_set)/1e3,y);
								
									Q = m[11]+m[2]-m[10]-m[3];
									Eex = Ein + Q - E_par - E_Nuc;
									
									Triton_Eex->Fill(Eex/1.e3);
									
									// ADVANCED
									// Calculating the offset
									
									Eex = 0;
									Ein = E_Nuc + E_par - Q - Eex;
									P_in = sqrt(2*m[11]*Ein);
									
									Triton_XvsdeltaEin1->Fill((Ein-Ein_set)/1e3,x);
									Triton_YvsdeltaEin1->Fill((Ein-Ein_set)/1e3,y);
									Triton_XvsdeltaEin1_det[(*FBarrel)[j_bar].GetID()]->Fill((Ein-Ein_set)/1e3,x);
									Triton_YvsdeltaEin1_det[(*FBarrel)[j_bar].GetID()]->Fill((Ein-Ein_set)/1e3,y);
		
//									fthetain->SetParameters(x_par,y_par,z_par,P_par,x_Nuc,y_Nuc,z_Nuc,P_Nuc,P_in);
//									x = fthetain->GetMinimumX(-10.,10);
									g->SetParameters(x_par,y_par,z_par,P_par,x_Nuc,y_Nuc,z_Nuc,P_Nuc,P_in,0,y);
									x = g->GetMinimumX(-10.,10.);
									h->SetParameters(x_par,y_par,x_Nuc,y_Nuc);
									y = h->Eval(x);
									fthetain->SetParameters(x_par,y_par,z_par,P_par,x_Nuc,y_Nuc,z_Nuc,P_Nuc,P_in,y);
									theta_in = atan(fthetain->Eval(x));
									g->SetParameters(x_par,y_par,z_par,P_par,x_Nuc,y_Nuc,z_Nuc,P_Nuc,P_in,x,y);
									z = g->GetMinimumX(-5.,5.);

									Triton_XY_2->Fill(x,y);
									Triton_Average_XY_2->Fill(fAveragePhix->Eval(x,y),fAveragePhiy->Eval(x,y));
									if (g->Eval(z)<1.) Triton_XY_3->Fill(x,y);
									Triton_Z->Fill(z);
									Triton_Ein1->Fill(Ein/1.e3);

									// Calculating the incoming momentum
									R_fun->SetParameters(x_Nuc,y_Nuc,z_Nuc);
									R_Nuc = R_fun->Eval(x,y,z);
									R_fun->SetParameters(x_par,y_par,z_par);
									R_par = R_fun->Eval(x,y,z);

									P_in_x = P_par*(x_par+x)/R_par+P_Nuc*(x_Nuc+x)/R_Nuc;	
									P_in_y = P_par*(y_par+y)/R_par+P_Nuc*(y_Nuc+y)/R_Nuc;	
									P_in_z = P_par*(z_par+z)/R_par+P_Nuc*(z_Nuc+z)/R_Nuc;	
							
									// Calculating the incoming angle
									if (P_in_x>=0&&P_in_y>=0) phi_in = atan(P_in_y/P_in_x);
									if (P_in_x>=0&&P_in_y<0) phi_in = 2*pi + atan(P_in_y/P_in_x);
									if (P_in_x<0&&P_in_y>=0) phi_in = pi + atan(P_in_y/P_in_x);
									if (P_in_x<0&&P_in_y<0) phi_in = pi + atan(P_in_y/P_in_x);

									Triton_Thetain_2->Fill(theta_in*1e3);
									Triton_Thetain_Phiin_2->Fill(theta_in*cos(phi_in)*1e3,theta_in*sin(phi_in)*1e3);
								}
							}
						}
					}
				}	
			}
		}
		// End forward barrel


		if(i%1000 == 0){
			cout<<setw(5)<<setiosflags(ios::fixed)<<setprecision(1)<<(100.*i)/nentries<<" % done\r"<<flush;
		}
	}

	// Writing histograms
	// Deuterons
	Deuteron_Average->Write("",TObject::kOverwrite);
	Deuteron_Average_XY->Write("",TObject::kOverwrite);
	Deuteron_XY->Write("",TObject::kOverwrite);
	Deuteron_Ein1->Write("",TObject::kOverwrite);
	Deuteron_Ein2->Write("",TObject::kOverwrite);
	Deuteron_EvsTheta->Write("",TObject::kOverwrite);
	Be11_EvsTheta->Write("",TObject::kOverwrite);
	Deuteron_dEvsEcor->Write("",TObject::kOverwrite);
	Deuteron_Thetain->Write("",TObject::kOverwrite);
	Deuteron_XvsThetain->Write("",TObject::kOverwrite);
	Deuteron_YvsThetain->Write("",TObject::kOverwrite);

	// Tritons
	Triton_Average->Write("",TObject::kOverwrite);
	Triton_Average_u2->Write("",TObject::kOverwrite);
	Triton_Average_XY->Write("",TObject::kOverwrite);
	Triton_Average_XY_2->Write("",TObject::kOverwrite);
	Triton_ZXpar->Write("",TObject::kOverwrite);
	Triton_ZYpar->Write("",TObject::kOverwrite);
	Triton_XY->Write("",TObject::kOverwrite);
	Triton_XY_2->Write("",TObject::kOverwrite);
	Triton_XY_3->Write("",TObject::kOverwrite);
	Triton_XY_Strip->Write("",TObject::kOverwrite);
	Triton_Z->Write("",TObject::kOverwrite);
	Triton_Ein1->Write("",TObject::kOverwrite);
	Triton_Ein2->Write("",TObject::kOverwrite);
	Triton_Eex->Write("",TObject::kOverwrite);
	Triton_dtvsEex->Write("",TObject::kOverwrite);
	Triton_EvsTheta->Write("",TObject::kOverwrite);
	Be10_EvsTheta->Write("",TObject::kOverwrite);
	
	
	For_dEvsEcor->Write("",TObject::kOverwrite);
	Triton_dEvsEcor->Write("",TObject::kOverwrite);
	Triton_Thetain->Write("",TObject::kOverwrite);
	Triton_Thetain_XY->Write("",TObject::kOverwrite);
	Triton_Thetain_2->Write("",TObject::kOverwrite);
	Triton_Thetain_Phiin_2->Write("",TObject::kOverwrite);
	Triton_XvsThetain->Write("",TObject::kOverwrite);
	Triton_YvsThetain->Write("",TObject::kOverwrite);
	Triton_XvsdeltaEin1->Write("",TObject::kOverwrite);
	Triton_YvsdeltaEin1->Write("",TObject::kOverwrite);
	Triton_XvsdeltaEin2->Write("",TObject::kOverwrite);
	Triton_YvsdeltaEin2->Write("",TObject::kOverwrite);
	
	for (Int_t idet=0;idet<4;idet++){
		Triton_XY_det[idet]->Write("",TObject::kOverwrite);
		Triton_Thetain_det[idet]->Write("",TObject::kOverwrite);
		Triton_XvsThetain_det[idet]->Write("",TObject::kOverwrite);
		Triton_YvsThetain_det[idet]->Write("",TObject::kOverwrite);
		Triton_XvsdeltaEin1_det[idet]->Write("",TObject::kOverwrite);
		Triton_YvsdeltaEin1_det[idet]->Write("",TObject::kOverwrite);
		Triton_XvsdeltaEin2_det[idet]->Write("",TObject::kOverwrite);
		Triton_YvsdeltaEin2_det[idet]->Write("",TObject::kOverwrite);
	}
	
	// Single histograms
	
	outfile->Close();
	delete tr;

	
	return 0;

}

